#ifndef CSUPBR_PBRMATH_H
#define CSUPBR_PBRMATH_H

#include "transform.h"

namespace csupbr {

    inline Vec3f reflect(const Vec3f &v, const Normal &n) {
        return v - 2 * dot(v, n) * n;
    }

    inline void reflect(const Vec3f &vi, Vec3f &vo, const Normal &n) {
        vo = reflect(vi, n);
    }

    inline bool refract(const Vec3f& in, Vec3f& out, const Normal& normal, float refRate) {
		Vec3f uvin = in.normalized();
		float dt = dot(uvin, normal);
		float discrimination = 1 - refRate * refRate * (1 - dt * dt);
		if (discrimination > 0) {
			out = refRate * (uvin - normal * dt) - normal * sqrt(discrimination);
			return true;
		}
		return false;
	}

    inline float schlick(float cosine, float refRate) {
		float r0 = (1.f - refRate) / (1.f + refRate);
		r0 *= r0;
		return r0 + (1.f - r0) * std::pow((1.f - cosine), 5.f);
	}

    [[deprecated]] inline Vec3f sphericalDirection(Float sinTheta, Float cosTheta, Float phi,
                                        const Vec3f &x, const Vec3f &y,
                                        const Vec3f &z) {
        return sinTheta * std::cos(phi) * x + sinTheta * std::sin(phi) * y +
                cosTheta * z;
    }

    [[deprecated]] inline Vec3f toWorld(const Vec3f &a, const Normal &N) {
        Vec3f B, C;
        if (std::abs(N.x) > std::abs(N.y)) {
            float invLen = 1.0f / std::sqrt(N.x * N.x + N.z * N.z);
            C = Vec3f(N.z * invLen, 0.0f, -N.x *invLen);
        } else {
            float invLen = 1.0f / std::sqrt(N.y * N.y + N.z * N.z);
            C = Vec3f(0.0f, N.z * invLen, -N.y *invLen);
        }
        B = cross(C, N);
        return a.x * B + a.y * C + a.z * N;
    }

    [[deprecated]] inline Vec3f random_in_unit_ball() {
		Vec3f re;
		do {
			re = 2.0f * Vec3f {
                (rand() % 1000) * .001f,
                (rand() % 1000) * .001f,
                (rand() % 1000) * .001f,
            } - Vec3f{1.f};
		} while (re.length() > 1.0f);
		return re;
	}

    inline Vec3f uniformSphereSampling(Point2f sample) {
        float theta = sample.x * Pi2;
        float phi = std::acos(1.0f - 2.0f * sample.y);
        return {
            std::sin(phi) * std::cos(theta),
            std::cos(phi),
            std::sin(phi) * std::sin(theta),
        };
    }

    // uniform Hemisphere Sampling
    inline Vec3f uniformHemisphereSampling(Point2f sample) {
        float phi = Pi2 * sample.y;
        float cosTheta = 1.f - sample.x;
        float sinTheta = std::sqrt(1.f - cosTheta * cosTheta);
        return {std::cos(phi) * sinTheta, std::sin(phi) * sinTheta, cosTheta};
    }

    // cosine Hemisphere Sampling
    inline Vec3f cosineHemisphereSampling(Point2f sample) {
        float phi = Pi2 * sample.y;
        float cosTheta = std::sqrt(1.f - sample.x);
        float sinTheta = std::sqrt(1.f - cosTheta * cosTheta);
        return {std::cos(phi) * sinTheta, std::sin(phi) * sinTheta, cosTheta};
    }

    inline Vec2f UniformSampleDisk(const Point2f &u) {
        Float r = std::sqrt(u[0]);
        Float theta = 2 * Pi * u[1];
        return Vec2f{r * std::cos(theta), r * std::sin(theta)};
    }

    inline Vec2f ConcentricDiskSample(const Point2f &rd) {
        Point2f uOffset = 2.f * rd - Point2f{1.f, 1.f};
        // Handle degeneracy at the origin
        if (uOffset.x == 0 && uOffset.y == 0)
            return {0, 0};
        // Apply concentric mapping to point
        Float theta, r;
        if (std::abs(uOffset.x) > std::abs(uOffset.y)) {
            r = uOffset.x;
            theta = Pi4 * (uOffset.y / uOffset.x);
        } else {
            r = uOffset.y;
            theta = Pi2 - Pi4 * (uOffset.x / uOffset.y);
        }
        return r * Vec2f{std::cos(theta), std::sin(theta)};
    }

    /**
     * 
     * \param n as z axis
     * \param s as x axis
     * \param t as y axis
     */
    inline void coordinateSystem(const Vec3f &n, Vec3f &s, Vec3f &t) {
        if (std::abs(n.x) > std::abs(n.y)) {
            Float invLen = 1.0f / std::sqrt(n.x * n.x + n.z * n.z);
            t = Vec3f(n.z * invLen, 0.0f, -n.x * invLen);
        } else {
            Float invLen = 1.0f / std::sqrt(n.y * n.y + n.z * n.z);
            t = Vec3f(0.0f, n.z * invLen, -n.y * invLen);
        }
        s = t.cross(n);
    }

    struct CoordinateSystem {
        explicit inline CoordinateSystem(const Normal &n = {.0f, .0f, 1.f}) : n(n) {
            coordinateSystem(n, s, t);
        }

        explicit inline CoordinateSystem(const Vec3f &s, const Vec3f &t, const Normal &n) : s(s), t(t), n(n) {}

        explicit inline CoordinateSystem(const Vec3f &x, const Vec3f &y, const Vec3f &z) : s(x), t(y), n(z) {}

        inline Vec3f toLocal(const Vec3f &v) {
            return {dot(v, s), dot(v, t), dot(v, n)};
        }

        inline Vec3f toWorld(const Vec3f &v) {
            return s * v.x + t * v.y + n * v.z;
        }

        Vec3f s, t;
        Normal n;
    }; 

	struct RNG {
		explicit RNG(uint64_t state = 0) : state(state + increment) {
			pcg32();
		}

		uint32_t uniformUint32() {
			return pcg32();
		}

		float uniformFloat() {
			return float(uniformUint32()) / std::numeric_limits<uint32_t>::max();
		}

    private:
        uint64_t state;
		static uint64_t const multiplier = 6364136223846793005u;
		static uint64_t const increment = 1442695040888963407u;

		static uint32_t rotr32(uint32_t x, unsigned r) {
			return x >> r | x << (-r & 31);
		}

		uint32_t pcg32() {
			uint64_t x = state;
			auto count = (unsigned)(x >> 59ULL);        // 59 = 64 - 5

			state = x * multiplier + increment;
			x ^= x >> 18ULL;                                // 18 = (64 - 27)/2
			return rotr32((uint32_t)(x >> 27ULL), count);    // 27 = 32 - 5
		}
		
	};

    inline Float FrDielectric(Float cosThetaI, Float etaI, Float etaT) {
        cosThetaI = clamp(cosThetaI, -1, 1);
        // <<Potentially swap indices of refraction>> 
        bool entering = cosThetaI > 0.f;
        if (!entering) {
            std::swap(etaI, etaT);
            cosThetaI = std::abs(cosThetaI);
        }

        // <<Compute cosThetaT using Snell’s law>> 
        Float sinThetaI = std::sqrt(std::max((Float)0,
                                                1 - cosThetaI * cosThetaI));
        Float sinThetaT = etaI / etaT * sinThetaI;
        //    <<Handle total internal reflection>> 
        Float cosThetaT = std::sqrt(std::max((Float)0,
                                                1 - sinThetaT * sinThetaT));

        Float Rparl = ((etaT * cosThetaI) - (etaI * cosThetaT)) /
                    ((etaT * cosThetaI) + (etaI * cosThetaT));
        Float Rperp = ((etaI * cosThetaI) - (etaT * cosThetaT)) /
                    ((etaI * cosThetaI) + (etaT * cosThetaT));
        return (Rparl * Rparl + Rperp * Rperp) / 2;
    }

}

#endif // !CSUPBR_PBRMATH_H